## ----measure-network-density----

load("replication_file_04_friendship_graph.RData")
load("replication_file_06_recruitment_graph.RData")
load('replication_file_07_geostat_grid.RData')

library(data.table)
library(parallel)
library(igraph)
library(sp)
library(sf)
library(rgdal)
library(dplyr)

region_2013.sp <- readOGR('.', 'Reg01012013_g_WGS84')
region_2013.sp <- 
  spTransform(region_2013.sp, "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
region_2013.sp$COD_REG <- as.numeric(as.character(region_2013.sp$COD_REG))

region_2013.sf <- st_as_sf(region_2013.sp)
region_south_2013.sf <- region_2013.sf %>% filter(COD_REG >= 13)

lonlat_node_df <- data.frame(lon = V(g_recruitment)$lon,
                             lat = V(g_recruitment)$lat)
lonlat_node.sp <- 
  SpatialPoints(lonlat_node_df, 
                proj4string=CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))
res <- over(lonlat_node.sp, region_2013.sp)

V(g_recruitment)$DEN_REG <- as.character(res$DEN_REG)
V(g_recruitment)$COD_REG <- res$COD_REG
V(g_recruitment)$north_south <- 
  ifelse(V(g_recruitment)$COD_REG < 13, "North", "South")

sequence <-
  seq(from = min(as.Date(V(g_recruitment)$joined_mu)), t
      o = as.Date(max(V(g_recruitment)$joined_mu)),
      by = 7)

geostat_grid.sp <- 
  as(geostat_grid_simp.sf,'Spatial')
geostat_grid.sp <- 
  spTransform(geostat_grid.sp, 
              '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0')

# Italy only
g_recruitment <- 
  g_recruitment - V(g_recruitment)[is.na(V(g_recruitment)$north_south)]
lonlat_node_df <- 
  data.frame(lon = V(g_recruitment)$lon,
             lat = V(g_recruitment)$lat)
lonlat_node.sp <- 
  SpatialPoints(lonlat_node_df, 
                proj4string=CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))

lonlat_north_node.df <- 
  data.frame(lon = V(g_recruitment)$lon[V(g_recruitment)$north_south == "North"],
             lat = V(g_recruitment)$lat[V(g_recruitment)$north_south == "North"])
lonlat_north_node.sp <- 
  SpatialPoints(lonlat_north_node.df, 
                proj4string=CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))

lonlat_south_node.df <- 
  data.frame(lon = V(g_recruitment)$lon[V(g_recruitment)$north_south == "South"],
             lat = V(g_recruitment)$lat[V(g_recruitment)$north_south == "South"])
lonlat_south_node.sp <- 
  SpatialPoints(as.matrix(lonlat_south_node.df), 
                proj4string=CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))

geostat_grid.sp <- 
  as(geostat_grid_simp.sf,'Spatial')
geostat_grid.sp <- 
  spTransform(geostat_grid.sp, 
              '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0')

res <- over(geostat_grid.sp, region_2013.sp)
geostat_grid.sp$COD_REG <- res$COD_REG
geostat_grid.sp$DEN_REG <- as.character(res$DEN_REG)

tot_pop <- sum(geostat_grid.sp$pop)
tot_north_pop <- sum(geostat_grid.sp$pop[geostat_grid.sp$COD_REG < 13], na.rm = T)
tot_south_pop <- sum(geostat_grid.sp$pop[geostat_grid.sp$COD_REG >= 13], na.rm = T)

## Tot
res <- over(lonlat_node.sp, geostat_grid.sp)
res <- as.data.frame(table(res$GRID_ID_simp))
geostat_grid.sp$n_nodes <- res$Freq[match(geostat_grid.sp$GRID_ID_simp, res$Var1)]
geostat_grid.sp$n_nodes[is.na(geostat_grid.sp$n_nodes)] <- 0

## North
res <- over(lonlat_north_node.sp, geostat_grid.sp)
res <- as.data.frame(table(res$GRID_ID_simp))
geostat_grid.sp$n_north_nodes <- res$Freq[match(geostat_grid.sp$GRID_ID_simp, res$Var1)]
geostat_grid.sp$n_north_nodes[is.na(geostat_grid.sp$n_north_nodes)] <- 0

## South
res <- over(lonlat_south_node.sp, geostat_grid.sp)
res <- as.data.frame(table(res$GRID_ID_simp))
geostat_grid.sp$n_south_nodes <- res$Freq[match(geostat_grid.sp$GRID_ID_simp, res$Var1)]
geostat_grid.sp$n_south_nodes[is.na(geostat_grid.sp$n_south_nodes)] <- 0

# Tot node
tot_max_act_cells <- sum(geostat_grid.sp$n_nodes > 0)
north_max_act_cells <- sum(geostat_grid.sp$n_north_nodes > 0)
south_max_act_cells <- sum(geostat_grid.sp$n_south_nodes > 0)

getPnt <- function(g) {
  require(igraph)
  lonlat_node.df <- 
    data.frame(lon = V(g)$lon,
               lat = V(g)$lat)
  lonlat_node.sp <- 
    SpatialPoints(lonlat_node.df,
                  proj4string=CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))
  return(lonlat_node.sp)
}

computeNetStats <- function(i) {
  require(igraph)
  require(sp)
  print(i)
  this_g <- 
    g_recruitment - V(g_recruitment)[as.Date(V(g_recruitment)$joined_mu) >= sequence[i]]
  
  this_g <- this_g -  V(this_g)[is.na(V(this_g)$north_south)]
  this_g_north <- this_g - V(this_g)[V(this_g)$north_south == 'South']
  this_g_south <- this_g - V(this_g)[V(this_g)$north_south == 'North']
  
  res <- over(getPnt(this_g), geostat_grid.sp)
  res <- unique(res[,c('GRID_ID_simp','pop')])
  this_act_cells <- nrow(res)
  this_act_pop <- sum(res$pop, na.rm = T)
  
  if (vcount(this_g_north) > 0) {
    res <- over(getPnt(this_g_north), geostat_grid.sp)
    res <- unique(res[,c('GRID_ID_simp','pop')])
    this_north_act_cells <- nrow(res)
    this_north_act_pop <- sum(res$pop, na.rm = T)
  } else {
    this_north_act_cells <- 0
    this_north_act_pop <- 0
  }
  
  res <- over(getPnt(this_g_south), geostat_grid.sp)
  res <- unique(res[,c('GRID_ID_simp','pop')])
  this_south_act_cells <- nrow(res)
  this_south_act_pop <- sum(res$pop, na.rm = T)
  
  density_cum_df <- 
    data.frame(time = sequence[i],
               density = edge_density(this_g),
               n = vcount(this_g),
               act_cells =  this_act_cells / tot_max_act_cells,
               act_pop = this_act_pop / tot_pop,
               scope = 'all')
  density_cum_df <- 
    rbind(density_cum_df, 
          data.frame(time = sequence[i],
                     density = edge_density(this_g_north),
                     n = vcount(this_g_north),
                     act_cells = this_north_act_cells / north_max_act_cells,
                     act_pop = this_north_act_pop / tot_north_pop,
                     scope = 'north'))
  density_cum_df <- 
    rbind(density_cum_df, 
          data.frame(time = sequence[i],
                     density = edge_density(this_g_south),
                     n = vcount(this_g_south),
                     act_cells = this_south_act_cells / south_max_act_cells,
                     act_pop = this_south_act_pop / tot_south_pop,
                     scope = 'south'))
  return(density_cum_df)
}

no_cores <- detectCores() - 1
cl <- makeCluster(no_cores, type="FORK")

density_cum_dt <- rbindlist(parLapply(cl, 2:length(sequence), computeNetStats))
stopCluster(cl)

save(density_cum_dt, file = "replication_file_09_network_density.RData")